clear all;
clc;

dXmin = 1000;
dXmax = 1000;
dYmin = 
dYmax = 
uiNumX = 
uiNumY = 
X = linspace(dXmin, dXmax, uiNumX);
Y = linspace(dYmin, dYmax, uiNumY); 

dZ = 0;
dA = 45; %x轴的正方向

%场源参数

dXc = (dXmin+dXmax)/2.0;
dYc = (dYmin+dYmax)/2.0;
dZc = 250;
dR = 100.0;

dK = 0.01;
dMr = 0.0;
dDr = 135.0; %度数
dIr = 45.0;

%背景场参数

dBO = 50000;
dI = 45.0;
dD = 0.0;

% 计算磁异常

DB1 = zeros(uiNumY, uiNumX);
%DB2 = zeros(uiNumY, uiNumX);

for i=1:uiNumY
	for j=1:uiNumX
		[DB1(i,j), DB2(i,j)] = sphere_DB(X(i), Y(i), dZ,...    % 观察点坐标
                      			dXc, dYc, dZc,... %球心坐标
                      			dR,...       %球体半径
                      			dK,...      %磁化率
                      			dMr, dIr, dDr,... %剩余磁化强度的大小， 倾角和偏角
                      			dA,...    %x轴正方向的方位角
                      			dBO, dI, dD);
	end
end

subplot(2,1,1);
contourf(X, Y, DB1, 20); colorbar; axis equal; set(gca, 'Ydir', 'reverse');
subplot(2,1,2);
contourf(X, Y, DB2, 20); colorbar; axis equal; set(gca, 'Ydir', 'reverse');
